Spatial Scripts

packrat::init("~/Desktop/Heloderma Spatial/Heloderma Spatial")

library(adehabitatHR) #for home range calculations
library(data.table) #manipulate S3 and S4 data tables

data.table 1.11.8  Latest news: r-datatable.com

Attaching package: ‘data.table’

The following objects are masked from ‘package:dplyr’:

    between, first, last

The following object is masked from ‘package:purrr’:

    transpose

The following object is masked from ‘package:raster’:

    shift

library(ggplot2) #for graphic output
library(ggfortify) #to allow ggplot2 to read spatial data
library(grid) #to add annotations to the output
# library(OpenStreetMap) #for obtaining raster images
library(pbapply) #needed for progress bar

package ‘pbapply’ was built under R version 3.5.2

library(plotly) #for interactive xy plot

Attaching package: ‘plotly’

The following objects are masked from ‘package:plyr’:

    arrange, mutate, rename, summarise

The following object is masked from ‘package:ggmap’:

    wind

The following object is masked from ‘package:ordinal’:

    slice

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:MASS’:

    select

The following object is masked from ‘package:raster’:

    select

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

library(rgdal) #for converting spatial data

package ‘rgdal’ was built under R version 3.5.2rgdal: version: 1.4-3, (SVN revision 828)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
 Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
 GDAL binary built with GEOS: FALSE 
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
 Linking to sp version: 1.3-1 

library(sp) #for converting spatial data
library(rgeos)

rgeos version: 0.4-2, (SVN revision 581)
 GEOS runtime version: 3.6.1-CAPI-1.10.1 
 Linking to sp version: 1.3-1 
 Polygon checking: TRUE 

# library(raster)
library(mapview)

Overall individual yearly home ranges for non/subsidized populations

Gila monster locations for all tracked lizards across Stone Canyon

All.Gilas <- read_csv("./GM_Final_Data.csv")

Missing column names filled in: 'X9' [9]Parsed with column specification:
cols(
  LIZARDNUMBER = col_character(),
  DATE = col_character(),
  YEAR = col_double(),
  SEASON = col_character(),
  TIME = col_character(),
  EASTING = col_double(),
  NORTHING = col_double(),
  `Distance Moved (m)` = col_double(),
  X9 = col_logical()
)

utm_points <- cbind(All.Gilas$EASTING, All.Gilas$NORTHING)
utm_locations <- SpatialPoints(utm_points, proj4string=CRS.SC)
proj_lat.lon <- as.data.frame(spTransform(utm_locations, CRS("+proj=longlat +datum=WGS84")))
colnames(proj_lat.lon) <- c("x","y")
## FORTIGY SPATIAL SPATIAL POINTS FOR PLOTTING:
proj_lat.lon <- fortify(proj_lat.lon, region = "Type")
myMap <- get_stamenmap(bbox = c(left = -111.009,
                                bottom = 32.459,
                                right = -110.969,
                                top = 32.474),
                       maptype = "terrain", 
                       crop = FALSE,
                       zoom = 15)

Source : http://tile.stamen.com/terrain/15/6279/13255.png
Source : http://tile.stamen.com/terrain/15/6280/13255.png
Source : http://tile.stamen.com/terrain/15/6281/13255.png
Source : http://tile.stamen.com/terrain/15/6282/13255.png
Source : http://tile.stamen.com/terrain/15/6283/13255.png
Source : http://tile.stamen.com/terrain/15/6279/13256.png
Source : http://tile.stamen.com/terrain/15/6280/13256.png
Source : http://tile.stamen.com/terrain/15/6281/13256.png
Source : http://tile.stamen.com/terrain/15/6282/13256.png
Source : http://tile.stamen.com/terrain/15/6283/13256.png
Source : http://tile.stamen.com/terrain/15/6279/13257.png
Source : http://tile.stamen.com/terrain/15/6280/13257.png
Source : http://tile.stamen.com/terrain/15/6281/13257.png
Source : http://tile.stamen.com/terrain/15/6282/13257.png
Source : http://tile.stamen.com/terrain/15/6283/13257.png

ggmap(myMap)+geom_point(data=proj_lat.lon, aes(x=x, y=y), size=0.3)

OVERALL YEARLY ANALYSES

Plot of 100% MCP HRs against number of relocations

year <- read_csv("GM_Consolidated_ByYear.csv")

Missing column names filled in: 'X13' [13]Parsed with column specification:
cols(
  Year = col_double(),
  Gila = col_character(),
  Sex = col_character(),
  Environment = col_character(),
  Home_Range_100mcp = col_double(),
  N100 = col_double(),
  Home_Range_95mcp = col_double(),
  N95 = col_double(),
  Home_Range_95kde = col_double(),
  N = col_double(),
  Home_Range_50kde = col_double(),
  N50 = col_double(),
  X13 = col_logical()
)

# quick plot
# Graph1<-ggplot(year,aes(x=N100,y=Home_Range_100mcp,group=Environment))+
Graph1<-ggplot(year,aes(x=N100,y=Home_Range_100mcp))+
  geom_point(aes(shape = factor(Environment)), size = 2)+
  geom_smooth(aes(linetype=Environment),colour="black", method="lm") +
  # scale_colour_manual(values=c(subsidized="cyan3",nonsubsidized="indian red1"))+
  # labs(title = "100% MCP Home Ranges")+
  xlab("Number of Relocations")+
  ylab("100% MCP Area (ha)")+
  # labs(caption = "Figure 3 | Non-Subsidized (Owl Head Buttes) vs. Subsidized (Stone Canyon) population 100% MCPs against number \n of fixes of the complete data set.")+
  theme(plot.caption = element_text(hjust = 0,lineheight = 0.9))
  # theme_bw()
Graph1<-Graph1+theme(axis.title=element_text(size = 14))
# legend at top-left, inside the plot
SCOH.hr.fig<-Graph1 + theme(legend.title = element_blank(),
               legend.text = element_text(size = 12),
               legend.justification=c(0,1),
               legend.position=c(0.05, 0.95),
               legend.background = element_blank(),
               legend.key = element_blank(),
               legend.box.background = element_rect(colour = "black")) +
               scale_shape_discrete(name  ="",
                          breaks=c("nonsubsidized", "subsidized"),
                          labels=c("Nonsubsidized", "Subsidized")) +
                            scale_linetype_discrete(name  ="",
                          breaks=c("nonsubsidized", "subsidized"),
                          labels=c("Nonsubsidized", "Subsidized"))
SCOH.hr.fig

# dir.create("outputs") # create a new folder to hold the output files
# ggsave("outputs/SC_OHB_plot.pdf")

Overall combined 100% MCP means averaged across sex

library(Rmisc)
Means <- summarySE(year, measurevar="Home_Range_100mcp",
                          groupvars=c("Environment"),na.rm = TRUE)
kable(Means, format = "pandoc", caption = 'Overall combined 100% MCP means averaged across sex')

Overall combined 100% MCP means averaged across sex

Environment	N	Home_Range_100mcp	sd	se	ci
nonsubsidized	26	33.44231	20.518658	4.0240400	8.287665
subsidized	53	10.40151	6.948743	0.9544832	1.915311

Overall combined 95% MCP means averaged across sex

Means.95mcp <- summarySE(year, measurevar="Home_Range_95mcp",
                          groupvars=c("Environment"),na.rm = TRUE)
Means.95mcp

Set projection for mapping

CRS.SC<-CRS("+proj=utm +zone=12 +ellps=WGS84 +units=m +no_defs")

Function for MCP analysis

Function of MCP polygons used for mapping

Function of KDE analysis

kde_analysis.href.plot <- function(filename, percentage){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$EASTING)
  y <- as.data.frame(data$NORTHING)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS.SC)
  xy <- SpatialPoints(data.proj@coords)
  kde<-kernelUD(xy, h="href", kern="bivnorm", grid=1000)
  ver <- getverticeshr(kde, percentage)
  area <- as.data.frame(round(ver$area,4))
  .rowNamesDF(area, make.names=TRUE) <- data$LIZARDNUMBER
  write.table(area,file="KDE_Hectares.csv",
              append=TRUE,sep=",", col.names=FALSE, row.names=TRUE)
  kde.points <- cbind((data.frame(data.proj@coords)),data$LIZARDNUMBER)
  colnames(kde.points) <- c("x","y","lizardnumber")
  kde.poly <- fortify(ver, region = "id")
  units <- grid.text(paste(round(ver$area,2)," ha"), x=0.9,  y=0.95,
                     gp=gpar(fontface=4, cex=0.9), draw = FALSE)
  kde.plot <- ggplot() +
    geom_polygon(data=kde.poly, aes(x=kde.poly$long, y=kde.poly$lat), alpha = 0.5) +
    geom_point(data=kde.points, aes(x=x, y=y)) + theme_bw() +
    labs(x="Easting (m)", y="Northing (m)", title=kde.points$lizardnumber) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) +
    annotation_custom(units)
  kde.plot
}

Function of KDE polygons for mapping

kde_analysis.href.polygon <- function(filename, percentage){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$EASTING)
  y <- as.data.frame(data$NORTHING)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS.SC)
  xy <- SpatialPoints(data.proj@coords)
  kde<-kernelUD(xy, h="href", kern="bivnorm", grid=1000)
  ver <- getverticeshr(kde, percentage)
  ver@proj4string<-CRS.SC
  area <- as.data.frame(round(ver$area,4))
  .rowNamesDF(area, make.names=TRUE) <- data$YEAR
  write.table(area,file="KDE_Hectares.csv",
              append=TRUE,sep=",", col.names=FALSE, row.names=TRUE)
  kde.points <- cbind((data.frame(data.proj@coords)),data$YEAR)
  colnames(kde.points) <- c("x","y","year")
  kde.poly <- fortify(ver, region = "id")
  units <- grid.text(paste(round(ver$area,2)," ha"), x=0.9,  y=0.95,
                     gp=gpar(fontface=4, cex=0.9), draw = FALSE)
  ver
}

Function of raster of UD

# kde_analysis.href.raster <- function(filename){
#   data <- read.csv(file = filename)
#   x <- as.data.frame(data$EASTING)
#   y <- as.data.frame(data$NORTHING)
#   xy <- c(x,y)
#   data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS.SC)
#   xy <- SpatialPoints(data.proj@coords)
#   kde<-kernelUD(xy, h="href", kern="bivnorm", grid=1000)
#   kde<-as(kde, "SpatialGridDataFrame")
#   kde@proj4string<- CRS.SC
#   kde
# }

Function of trajectory analysis and distance over time

traj_analysis <- function(filename){
  relocs_data <- read.csv(file = filename)
  relocs <- as.ltraj(cbind(relocs_data$EASTING, relocs_data$NORTHING),id=relocs_data$LIZARDNUMBER, typeII = FALSE, date=NULL)
  relocs.df <- ld(relocs)
  relocs_dist <- as.data.frame(sum(sapply(relocs.df$dist, sum, na.rm=TRUE)))
  colnames(relocs_dist) <- "Total Distance"
  name <- relocs.df$id[1]
  row.names(relocs_dist) <- name
  relocs_units <- grid.text(paste(round(relocs_dist,2),"m"), x=0.9, y=0.9,
                            gp=gpar(fontface=3, col="black", cex=0.9), draw = FALSE)
  reloc.plot <- ggplot() + theme_classic() + geom_path(data=relocs.df, aes(x=x,y=y), linetype = "dashed", colour = "red",
                                                       arrow = arrow(length=unit(.5,"cm"), angle = 20, ends="last", type = "closed")) +
    geom_point(data=relocs.df, aes(x=x, y=y)) + geom_point(data=relocs.df, aes(x=x[1],
                                                                               y=y[1]), size = 3, color = "darkgreen", pch=0) +
    labs(x="Easting (m)", y="Northing (m)", title=relocs.df$id[1]) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) +
    annotation_custom(relocs_units)
  reloc.plot
}

Function of distance of time

dist_analysis <- function(filename){
  relocs_data <- read.csv(file = filename)
  relocs <- as.ltraj(cbind(relocs_data$EASTING, relocs_data$NORTHING),id=relocs_data$LIZARDNUMBER, typeII = FALSE, date=NULL)
  relocs.df <- ld(relocs)
  relocs_dist <- as.data.frame(sum(sapply(relocs.df$dist, sum, na.rm=TRUE)))
  colnames(relocs_dist) <- "Total Distance"
  name <- relocs.df$id[1]
  row.names(relocs_dist) <- name
  write.table(relocs_dist,file="reloc_dist.csv",
              append=TRUE,sep=",", col.names=FALSE, row.names=TRUE)
  dist.plot
}

Map of yearly HR shifts of a subset of Gila Monsters. Includes running MCP polygons, Fortify mcp polygons for ggplot2 by YEAR

M215_mcp.11<-mcp_analysis.POLY("./M215/2011 .csv", percentage= 100)
M215_mcp.12<-mcp_analysis.POLY("./M215/2012 .csv", percentage= 100)
F104_mcp.08<-mcp_analysis.POLY("./F104/2008 .csv", percentage= 100)
F104_mcp.09<-mcp_analysis.POLY("./F104/2009 .csv", percentage= 100)
F114_mcp.08<-mcp_analysis.POLY("./F114/2008 .csv", percentage= 100)
F114_mcp.09<-mcp_analysis.POLY("./F114/2009 .csv", percentage= 100)
F114_mcp.10<-mcp_analysis.POLY("./F114/2010 .csv", percentage= 100)
F114_mcp.11<-mcp_analysis.POLY("./F114/2011 .csv", percentage= 100)
F114_mcp.12<-mcp_analysis.POLY("./F114/2012 .csv", percentage= 100)
F137_mcp.09<-mcp_analysis.POLY("./F137/2009 .csv", percentage= 100)
F137_mcp.10<-mcp_analysis.POLY("./F137/2010 .csv", percentage= 100)
F137_mcp.11<-mcp_analysis.POLY("./F137/2011 .csv", percentage= 100)
F147_mcp.09<-mcp_analysis.POLY("./F147/2009 .csv", percentage= 100)
F147_mcp.10<-mcp_analysis.POLY("./F147/2010 .csv", percentage= 100)
F147_mcp.11<-mcp_analysis.POLY("./F147/2011 .csv", percentage= 100)
F147_mcp.12<-mcp_analysis.POLY("./F147/2012 .csv", percentage= 100)
F36_mcp.08<-mcp_analysis.POLY("./F36/2008 .csv", percentage= 100)
F36_mcp.09<-mcp_analysis.POLY("./F36/2009 .csv", percentage= 100)
F36_mcp.10<-mcp_analysis.POLY("./F36/2010 .csv", percentage= 100)
F36_mcp.11<-mcp_analysis.POLY("./F36/2011 .csv", percentage= 100)
F36_mcp.12<-mcp_analysis.POLY("./F36/2012 .csv", percentage= 100)
F66_mcp.08<-mcp_analysis.POLY("./F66/2008 .csv", percentage= 100)
F66_mcp.09<-mcp_analysis.POLY("./F66/2009 .csv", percentage= 100)
F66_mcp.10<-mcp_analysis.POLY("./F66/2010 .csv", percentage= 100)
M119_mcp.08<-mcp_analysis.POLY("./M119/2008 .csv", percentage= 100)
M119_mcp.09<-mcp_analysis.POLY("./M119/2009 .csv", percentage= 100)
M119_mcp.10<-mcp_analysis.POLY("./M119/2010 .csv", percentage= 100)
M112_mcp.07<-mcp_analysis.POLY("./M112/2007 .csv", percentage= 100)
M112_mcp.09<-mcp_analysis.POLY("./M112/2009 .csv", percentage= 100)
M112_mcp.10<-mcp_analysis.POLY("./M112/2010 .csv", percentage= 100)
M69_mcp.09<-mcp_analysis.POLY("./M69/2009 .csv", percentage= 100)
M69_mcp.10<-mcp_analysis.POLY("./M69/2010 .csv", percentage= 100)

## Fortify mcp polygons for ggplot2 *YEAR*:

F104_mcp.08T <- fortify(F104_mcp.08, region = "id")
F104_mcp.09T <- fortify(F104_mcp.09, region = "id")
F114_mcp.08T <- fortify(F114_mcp.08, region = "id")
F114_mcp.09T <- fortify(F114_mcp.09, region = "id")
F114_mcp.10T <- fortify(F114_mcp.10, region = "id")
F114_mcp.11T <- fortify(F114_mcp.11, region = "id")
F114_mcp.12T <- fortify(F114_mcp.12, region = "id")
F137_mcp.09T <- fortify(F137_mcp.09, region = "id")
F137_mcp.10T <- fortify(F137_mcp.10, region = "id")
F137_mcp.11T <- fortify(F137_mcp.11, region = "id")
F147_mcp.09T <- fortify(F147_mcp.09, region = "id")
F147_mcp.10T <- fortify(F147_mcp.10, region = "id")
F147_mcp.11T <- fortify(F147_mcp.11, region = "id")
F147_mcp.12T <- fortify(F147_mcp.12, region = "id")
F36_mcp.08T <- fortify(F36_mcp.08, region = "id")
F36_mcp.09T <- fortify(F36_mcp.09, region = "id")
F36_mcp.10T <- fortify(F36_mcp.10, region = "id")
F36_mcp.11T <- fortify(F36_mcp.11, region = "id")
F36_mcp.12T <- fortify(F36_mcp.12, region = "id")
F66_mcp.08T <- fortify(F66_mcp.08, region = "id")
F66_mcp.09T <- fortify(F66_mcp.09, region = "id")
F66_mcp.10T <- fortify(F66_mcp.10, region = "id")
M119_mcp.08T <- fortify(M119_mcp.08, region = "id")
M119_mcp.09T <- fortify(M119_mcp.09, region = "id")
M119_mcp.10T <- fortify(M119_mcp.10, region = "id")
M112_mcp.07T <- fortify(M112_mcp.07, region = "id")
M112_mcp.09T <- fortify(M112_mcp.09, region = "id")
M112_mcp.10T <- fortify(M112_mcp.10, region = "id")
M69_mcp.09T <- fortify(M69_mcp.09, region = "id")
M69_mcp.10T <- fortify(M69_mcp.10, region = "id")
M215_mcp.11T <- fortify(M215_mcp.11, region = "id")
M215_mcp.12T <- fortify(M215_mcp.12, region = "id")


mcp.shift.TEST4 <- ggplot() +
  # geom_polygon(data=F104_mcp.08T, aes(x=F104_mcp.08T$long, y=F104_mcp.08T$lat),
  #              alpha=0.1,colour="black",linetype=2) +
  # geom_polygon(data=F104_mcp.09T, aes(x=F104_mcp.09T$long, y=F104_mcp.09T$lat),
  #              alpha=0.1,colour="black",linetype=2) +
  geom_polygon(data=F114_mcp.08T, aes(x=F114_mcp.08T$long, y=F114_mcp.08T$lat),
               alpha=0.1,colour="black",linetype=3) +
  geom_polygon(data=F114_mcp.09T, aes(x=F114_mcp.09T$long, y=F114_mcp.09T$lat),
               alpha=0.1,colour="black",linetype=3) +
  geom_polygon(data=F114_mcp.10T, aes(x=F114_mcp.10T$long, y=F114_mcp.10T$lat),
               alpha=0.1,colour="black",linetype=3) +
  geom_polygon(data=F114_mcp.11T, aes(x=F114_mcp.11T$long, y=F114_mcp.11T$lat),
               alpha=0.1,colour="black",linetype=3) +
  geom_polygon(data=F114_mcp.12T, aes(x=F114_mcp.12T$long, y=F114_mcp.12T$lat),
               alpha=0.1,colour="black",linetype=3) +
  geom_polygon(data=F137_mcp.09T, aes(x=F137_mcp.09T$long, y=F137_mcp.09T$lat),
               alpha=0.1,colour="brown",linetype=4) +
  geom_polygon(data=F137_mcp.10T, aes(x=F137_mcp.10T$long, y=F137_mcp.10T$lat),
               alpha=0.1,colour="brown",linetype=4) +
  geom_polygon(data=F137_mcp.11T, aes(x=F137_mcp.11T$long, y=F137_mcp.11T$lat),
               alpha=0.1,colour="brown",linetype=4) +
  geom_polygon(data=F147_mcp.09T, aes(x=F147_mcp.09T$long, y=F147_mcp.09T$lat),
               alpha=0.1,colour="red",linetype=1) +
  geom_polygon(data=F147_mcp.10T, aes(x=F147_mcp.10T$long, y=F147_mcp.10T$lat),
               alpha=0.1,colour="red",linetype=1) +
  geom_polygon(data=F147_mcp.11T, aes(x=F147_mcp.11T$long, y=F147_mcp.11T$lat),
               alpha=0.1,colour="red",linetype=1) +
  geom_polygon(data=F147_mcp.12T, aes(x=F147_mcp.12T$long, y=F147_mcp.12T$lat),
               alpha=0.1,colour="red",linetype=1) +
  # geom_polygon(data=F36_mcp.08T, aes(x=F36_mcp.08T$long, y=F36_mcp.08T$lat),
  #              alpha=0.1,colour="black",linetype=6) +
  # geom_polygon(data=F36_mcp.09T, aes(x=F36_mcp.09T$long, y=F36_mcp.09T$lat),
  #              alpha=0.1,colour="black",linetype=6) +
  # geom_polygon(data=F36_mcp.10T, aes(x=F36_mcp.10T$long, y=F36_mcp.10T$lat),
  #              alpha=0.1,colour="black",linetype=6) +
  # geom_polygon(data=F36_mcp.11T, aes(x=F36_mcp.11T$long, y=F36_mcp.11T$lat),
  #              alpha=0.1,colour="black",linetype=6) +
  # geom_polygon(data=F36_mcp.12T, aes(x=F36_mcp.12T$long, y=F36_mcp.12T$lat),
  #              alpha=0.1,colour="black",linetype=6) +
  geom_polygon(data=F66_mcp.08T, aes(x=F66_mcp.08T$long, y=F66_mcp.08T$lat),
               alpha=0.1,colour="green",linetype=5) +
  geom_polygon(data=F66_mcp.09T, aes(x=F66_mcp.09T$long, y=F66_mcp.09T$lat),
               alpha=0.1,colour="green",linetype=5) +
  geom_polygon(data=F66_mcp.10T, aes(x=F66_mcp.10T$long, y=F66_mcp.10T$lat),
               alpha=0.1,colour="green",linetype=5) +
  geom_polygon(data=M119_mcp.08T, aes(x=M119_mcp.08T$long, y=M119_mcp.08T$lat),
               alpha=0.1,colour="blue",linetype=6) +
  geom_polygon(data=M119_mcp.09T, aes(x=M119_mcp.09T$long, y=M119_mcp.09T$lat),
               alpha=0.1,colour="blue",linetype=6) +
  geom_polygon(data=M119_mcp.10T, aes(x=M119_mcp.10T$long, y=M119_mcp.10T$lat),
               alpha=0.1,colour="blue",linetype=6) +
  geom_polygon(data=M112_mcp.07T, aes(x=M112_mcp.07T$long, y=M112_mcp.07T$lat),
               alpha=0.1,colour="purple",linetype=2) +
  geom_polygon(data=M112_mcp.09T, aes(x=M112_mcp.09T$long, y=M112_mcp.09T$lat),
               alpha=0.1,colour="purple",linetype=2) +
  geom_polygon(data=M112_mcp.10T, aes(x=M112_mcp.10T$long, y=M112_mcp.10T$lat),
               alpha=0.1,colour="purple",linetype=2) +
  # geom_polygon(data=M69_mcp.09T, aes(x=M69_mcp.09T$long, y=M69_mcp.09T$lat),
  #              alpha=0.1,colour="black") +
  # geom_polygon(data=M69_mcp.10T, aes(x=M69_mcp.10T$long, y=M69_mcp.10T$lat),
  #              alpha=0.1,colour="black") +
  # geom_polygon(data=M215_mcp.11T, aes(x=M215_mcp.11T$long, y=M215_mcp.11T$lat),
  #              alpha=0.1,colour="black") +
  # geom_polygon(data=M215_mcp.12T, aes(x=M215_mcp.12T$long, y=M215_mcp.12T$lat),
  #              alpha=0.1,colour="black") +
  theme_bw() +labs(x="Easting (m)", y="Northing (m)") +
  labs(caption = "Figure 5  |  SC Yearly home range shifts of 8 lizards, both males and females. Home range shifts appear to be \n relativley stable over study years.")+
  theme(plot.caption = element_text(hjust = 0,lineheight = 0.9))
  # theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5))
## within each geom_polygon line?:
## aes(colour="red"or"M112_mcp.09T")...+scale_color_manual(name="",breaks=c("","",...""))+
## values=c(""="",...)

mcp.shift.TEST4

Raw group 100% MCP home range means of Stone Canyon and Owl Head Buttes. Grouped by environment and sex

library(Rmisc)
YR_GRP_Means <- summarySE(year, measurevar="Home_Range_100mcp",
                          groupvars=c("Environment","Sex"),na.rm = TRUE)

kable(YR_GRP_Means, format = "pandoc", 
      caption = 'Table 1 | Raw group 100% MCP home range means of Stone Canyon and Owl Head Buttes. Grouped by environment and sex.')

Raw group 95% MCP home range means of Stone Canyon and Owl Head Buttes. Grouped by environment and sex

YR_GRP_Means95 <- summarySE(year, measurevar="Home_Range_95mcp",
                            groupvars=c("Environment","Sex"),na.rm = TRUE)

kable(YR_GRP_Means95, format = "pandoc", caption = 'Table 2 | Raw group 95% MCP home range means of raw data of Stone Canyon and Owl Head Buttes. Grouped by environment and sex.')

RM-ANOVA for 100% MCP analyses between the subsidized and non-subsidized

# Get p-values from mixed model F values:
library(lme4)
library(readr)
year <- read_csv("GM_Consolidated_ByYear.csv")

RMmod.year<-lmer(Home_Range_100mcp~Environment+Year+Sex+N100+Environment*Sex+
                   (1|Gila),data = year)
summary(RMmod.year)

ANOVA table for 100% MCPs between the two populations

anova(RMmod.year)

RM-ANOVA for 95% MCP analyses between the subsidized and non-subsidized

RMmod.year95<-lmer(Home_Range_95mcp~Environment+Year+Sex+N100+Environment*Sex+
                   (1|Gila),data = year)
summary(RMmod.year95)

ANOVA table for 95% MCPs between the two populations

anova(RMmod.year95)

Raw group means AND adjusted EMMs of Yearly Overall 100%MCP

Directional means of home range (100% MCP) after being adjusted for year, sex and sample size

kable(ref_dfRM_sex, format = "pandoc", caption = 'Table | Directional means of home range (100% MCP) after being adjusted for year, sex and sample size.')

RM-ANOVA of 95% KDEs for the subsidized population

RM.KDEmod.year<-lmer(Home_Range_95kde~Year+Sex+N+(1|Gila),data = sub)

summary(RM.KDEmod.year)

ANOVA Table for 95%KDE

anova(RM.KDEmod.year)

RM-ANOVA of 50% KDEs for the subsidized

RM.KDE.50.mod.year<-lmer(Home_Range_50kde~Year+Sex+N+(1|Gila),data = sub)

summary(RM.KDE.50.mod.year)

ANOVA Talbe of 50% KDE for the subsidized

anova(RM.KDE.50.mod.year)

TABLE. Raw Group 50% KDE home range means male and female home ranges at Stone Canyon

YR_GRP_Means.50KDE <- summarySE(sub, measurevar="Home_Range_50kde",
                            groupvars=c("Sex"),na.rm = TRUE)

kable(YR_GRP_Means.50KDE, format = "pandoc", caption = 'Table 5 | Raw Group 50% KDE home range means male and female home ranges at Stone Canyon.')

Raw group means AND adjusted EMMs of Yearly Overall 95% KDEs between non/subsidized populations

Collective grid of 100% MCP and 95% KDE of both sites from above

ggarrange(Raw.YearHR, yr.mean.adj, Raw.kde, kde.mean.adj, labels = c("A", "B", "C","D"),
          ncol = 2, nrow = 2)

43.4 male 42.9 female Yearly overall means of 95% KDEs grouped by site and sex

YR_Means.95KDEall <- summarySE(year, measurevar="Home_Range_95kde",
                            groupvars=c("Environment","Sex"),na.rm = TRUE)
 
kable(YR_Means.95KDEall, format = "pandoc", caption = 'Table | Raw Group 95% KDE home range means male and female home ranges at non/subsidized.')

Table | Raw Group 95% KDE home range means male and female home ranges at non/subsidized.

Environment	Sex	N	Home_Range_95kde	sd	se	ci
nonsubsidized	female	5	36.80000	9.603905	4.294997	11.924824
nonsubsidized	male	6	69.40000	27.763789	11.334520	29.136310
subsidized	female	37	22.98892	11.046272	1.815996	3.683010
subsidized	male	13	35.00308	12.057546	3.344161	7.286302

Pairwise Comparisons, between sexes by environment, and between environments averaged across sex

RMmod.year.Em<-lmer(Home_Range_100mcp~Environment+Year+Sex+N100+Environment*Sex+
                      (1|Gila),data = year)

# RMmod.year.Em95<-lmer(Home_Range_95mcp~Environment+Year+Sex+N95+Environment*Sex+
#                       (1|Gila),data = year)

# Sex.emm.oa <- emmeans(RMmod.year.Em, c("Environment","Sex"))
# pairs(Sex.emm.oa)

emm_s.t2 <- emmeans(RMmod.year.Em, pairwise ~ Sex | Environment)
emm_s.t2
# emm_s.e1 <- emmeans(RMmod.year.Em, pairwise ~ Environment)
# emm_s.e1

Graphical Comparisons of Sex Within Each Environment:

plot(emm_s.t2, comparisons = TRUE, xlab = "Least Square Mean (ha)", ylab = "Environment")

Pairwise by sex between enviornements 100% MCP, and 95% KDEs

emm_s.t3 <- emmeans(RMmod.year.Em, pairwise ~ Environment | Sex)
emm_s.t3

# emm_s.t95 <- emmeans(RMmod.year.Em95, pairwise ~ Environment | Sex)
# emm_s.t95

Graphical Comparisons of Sex between the two populations:

plot(emm_s.t3, comparisons = TRUE, xlab = "Least Square Mean (ha)", ylab = "Environment")

Ineractive map of MCPs at Stone Canyon

M67_MCP<-mcp_analysis.POLY('./M67/M67 .csv', percentage= 100)
M69_MCP<-mcp_analysis.POLY('./M69/M69 .csv', percentage= 100)
M255_MCP<-mcp_analysis.POLY('./M255/M255 .csv', percentage= 100)
M215_MCP<-mcp_analysis.POLY('./M215/M215 .csv', percentage= 100)
M14_MCP<-mcp_analysis.POLY('./M14/M14 .csv', percentage= 100)
M119_MCP<-mcp_analysis.POLY('./M119/M119 .csv', percentage= 100)
M112_MCP<-mcp_analysis.POLY('./M112/M112 .csv', percentage= 100)

F66_MCP<-mcp_analysis.POLY('./F66/F66 .csv', percentage= 100)
F36_MCP<-mcp_analysis.POLY('./F36/F36 .csv', percentage= 100)
F252_MCP<-mcp_analysis.POLY('./F252/F252 .csv', percentage= 100)
F214_MCP<-mcp_analysis.POLY('./F214/F214 .csv', percentage= 100)
F200_MCP<-mcp_analysis.POLY('./F200/F200 .csv', percentage= 100)
F147_MCP<-mcp_analysis.POLY('./F147/F147 .csv', percentage= 100)
F146_MCP<-mcp_analysis.POLY('./F146/F146 .csv', percentage= 100)
F137_MCP<-mcp_analysis.POLY('./F137/F137 .csv', percentage= 100)
F135_MCP<-mcp_analysis.POLY('./F135/F135 .csv', percentage= 100)
F114_MCP<-mcp_analysis.POLY('./F114/F114 .csv', percentage= 100)
F104_MCP<-mcp_analysis.POLY('./F104/F104 .csv', percentage= 100)

Male.MCP <- rbind(M67_MCP,M69_MCP,M255_MCP,M215_MCP,M14_MCP,M119_MCP,M112_MCP)
Female.MCP <- rbind(F66_MCP,F36_MCP,F252_MCP,F214_MCP,F200_MCP,F147_MCP,F146_MCP,F137_MCP,
                    F135_MCP,F114_MCP,F104_MCP)

mapviewOptions(basemaps = c("OpenStreetMap","Esri.WorldImagery","OpenTopoMap"),
               na.color = "magenta",
               layers.control.pos = "topleft")

mapview(Male.MCP, legend=F, zcol="id", col.regions = c("blue"), alpha.regions=0.3) + 
  mapview(Female.MCP, legend=F, zcol = "id", col.regions = c("red"), alpha.regions=0.3)

Create stagnant stamen map of MCPs at Stone Canyon

Error: Don't know how to add north2(SC_stamen_map, x = 0.89, y = 0.85, scale = 0.1, symbol = 16) to a plot

Interactive map of KDEs at Stone Canyon

M67_KDE<-kde_analysis.href.polygon('./M67/M67 .csv', percentage= 95)
M69_KDE<-kde_analysis.href.polygon('./M69/M69 .csv', percentage= 95)
M255_KDE<-kde_analysis.href.polygon('./M255/M255 .csv', percentage= 95)
M215_KDE<-kde_analysis.href.polygon('./M215/M215 .csv', percentage= 95)
M14_KDE<-kde_analysis.href.polygon('./M14/M14 .csv', percentage= 95)
M119_KDE<-kde_analysis.href.polygon('./M119/M119 .csv', percentage= 95)
M112_KDE<-kde_analysis.href.polygon('./M112/M112 .csv', percentage= 95)

F66_KDE<-kde_analysis.href.polygon('./F66/F66 .csv', percentage= 95)
F36_KDE<-kde_analysis.href.polygon('./F36/F36 .csv', percentage= 95)
F252_KDE<-kde_analysis.href.polygon('./F252/F252 .csv', percentage= 95)
F214_KDE<-kde_analysis.href.polygon('./F214/F214 .csv', percentage= 95)
F200_KDE<-kde_analysis.href.polygon('./F200/F200 .csv', percentage= 95)
F147_KDE<-kde_analysis.href.polygon('./F147/F147 .csv', percentage= 95)
F146_KDE<-kde_analysis.href.polygon('./F146/F146 .csv', percentage= 95)
F137_KDE<-kde_analysis.href.polygon('./F137/F137 .csv', percentage= 95)
F135_KDE<-kde_analysis.href.polygon('./F135/F135 .csv', percentage= 95)
F114_KDE<-kde_analysis.href.polygon('./F114/F114 .csv', percentage= 95)
F104_KDE<-kde_analysis.href.polygon('./F104/F104 .csv', percentage= 95)

Male.KDE <- rbind(M67_KDE,M69_KDE,M255_KDE,M215_KDE,M14_KDE,M119_KDE,M112_KDE)
Female.KDE <- rbind(F66_KDE,F36_KDE,F252_KDE,F214_KDE,F200_KDE,F147_KDE,F146_KDE,F137_KDE,
                    F135_KDE,F114_KDE,F104_KDE)

mapviewOptions(basemaps = c("OpenStreetMap","Esri.WorldImagery","OpenTopoMap"),
               na.color = "magenta",
               layers.control.pos = "topleft")

mapview(Male.KDE, legend=F, zcol="id", col.regions = c("blue"), alpha.regions=0.3) + 
  mapview(Female.KDE, legend=F, zcol = "id", col.regions = c("red"), alpha.regions=0.3)

TABLE

Table | Subsidized and non-subsidized directional means of KDE home ranges after being adjusted for year, sex and sample size.

Environment	Sex	lsmean	SE	df	lower.CL	upper.CL
nonsubsidized	female	42.44837	9.346459	53.01148	23.70185	61.19490
subsidized	female	20.93992	4.083131	22.92042	12.49169	29.38814
nonsubsidized	male	80.81873	9.084540	52.50968	62.59348	99.04398
subsidized	male	35.27439	5.501543	29.24014	24.02648	46.52229

SEASONAL ANALYSES

Map of seasonal fluctions of home ranges

## Create MCP polygons by SEASON:
M215_mcp.EM<-mcp_analysis.POLY("./M215/Emergence .csv", percentage= 100)
M215_mcp.DRY<-mcp_analysis.POLY("./M215/Dry .csv", percentage= 100)
M215_mcp.MON<-mcp_analysis.POLY("./M215/Monsoon .csv", percentage= 100)

M112_mcp.DRY<-mcp_analysis.POLY("./M112/Dry .csv", percentage= 100)
M112_mcp.MON<-mcp_analysis.POLY("./M112/Monsoon .csv", percentage= 100)
M112_mcp.PM<-mcp_analysis.POLY("./M112/Post_Monsoon .csv", percentage= 100)

M119_mcp.DRY<-mcp_analysis.POLY("./M119/Dry .csv", percentage= 100)
M119_mcp.MON<-mcp_analysis.POLY("./M119/Monsoon .csv", percentage= 100)
M119_mcp.PM<-mcp_analysis.POLY("./M119/Post_Monsoon .csv", percentage= 100)

F114_mcp.EM<-mcp_analysis.POLY("./F114/Emergence .csv", percentage= 100)
F114_mcp.DRY<-mcp_analysis.POLY("./F114/Dry .csv", percentage= 100)
F114_mcp.MON<-mcp_analysis.POLY("./F114/Monsoon .csv", percentage= 100)
F114_mcp.PM<-mcp_analysis.POLY("./F114/Post_Monsoon .csv", percentage= 100)

F137_mcp.EM<-mcp_analysis.POLY("./F137/Emergence .csv", percentage= 100)
F137_mcp.DRY<-mcp_analysis.POLY("./F137/Dry .csv", percentage= 100)
F137_mcp.MON<-mcp_analysis.POLY("./F137/Monsoon .csv", percentage= 100)
F137_mcp.PM<-mcp_analysis.POLY("./F137/Post_Monsoon .csv", percentage= 100)

F147_mcp.EM<-mcp_analysis.POLY("./F147/Emergence .csv", percentage= 100)
F147_mcp.DRY<-mcp_analysis.POLY("./F147/Dry .csv", percentage= 100)
F147_mcp.MON<-mcp_analysis.POLY("./F147/Monsoon .csv", percentage= 100)
F147_mcp.PM<-mcp_analysis.POLY("./F147/Post_Monsoon .csv", percentage= 100)

F252_mcp.EM<-mcp_analysis.POLY("./F252/Emergence .csv", percentage= 100)
F252_mcp.DRY<-mcp_analysis.POLY("./F252/Dry .csv", percentage= 100)
F252_mcp.MON<-mcp_analysis.POLY("./F252/Monsoon .csv", percentage= 100)
F252_mcp.PM<-mcp_analysis.POLY("./F252/Post_Monsoon .csv", percentage= 100)

F36_mcp.EM<-mcp_analysis.POLY("./F36/Emergence .csv", percentage= 100)
F36_mcp.DRY<-mcp_analysis.POLY("./F36/Dry .csv", percentage= 100)
F36_mcp.MON<-mcp_analysis.POLY("./F36/Monsoon .csv", percentage= 100)
F36_mcp.PM<-mcp_analysis.POLY("./F36/Post_Monsoon .csv", percentage= 100)

F66_mcp.EM<-mcp_analysis.POLY("./F66/Emergence .csv", percentage= 100)
F66_mcp.DRY<-mcp_analysis.POLY("./F66/Dry .csv", percentage= 100)
F66_mcp.MON<-mcp_analysis.POLY("./F66/Monsoon .csv", percentage= 100)
F66_mcp.PM<-mcp_analysis.POLY("./F66/Post_Monsoon .csv", percentage= 100)

## Fortify mcp polygons for ggplot2 *SEASON*:
M215_mcp.EMT <- fortify(M215_mcp.EM, region = "id")
M215_mcp.DRYT <- fortify(M215_mcp.DRY, region = "id")
M215_mcp.MONT <- fortify(M215_mcp.MON, region = "id")

M112_mcp.DRYT <- fortify(M112_mcp.DRY, region = "id")
M112_mcp.MONT <- fortify(M112_mcp.MON, region = "id")
M112_mcp.PMT <- fortify(M112_mcp.PM, region = "id")

M119_mcp.DRYT <- fortify(M119_mcp.DRY, region = "id")
M119_mcp.MONT <- fortify(M119_mcp.MON, region = "id")
M119_mcp.PMT <- fortify(M119_mcp.PM, region = "id")

F114_mcp.EMT <- fortify(F114_mcp.EM, region = "id")
F114_mcp.DRYT <- fortify(F114_mcp.DRY, region = "id")
F114_mcp.MONT <- fortify(F114_mcp.MON, region = "id")
F114_mcp.PMT <- fortify(F114_mcp.PM, region = "id")

F137_mcp.EMT <- fortify(F137_mcp.EM, region = "id")
F137_mcp.DRYT <- fortify(F137_mcp.DRY, region = "id")
F137_mcp.MONT <- fortify(F137_mcp.MON, region = "id")
F137_mcp.PMT <- fortify(F137_mcp.PM, region = "id")

F147_mcp.EMT <- fortify(F147_mcp.EM, region = "id")
F147_mcp.DRYT <- fortify(F147_mcp.DRY, region = "id")
F147_mcp.MONT <- fortify(F147_mcp.MON, region = "id")
F147_mcp.PMT <- fortify(F147_mcp.PM, region = "id")

F252_mcp.EMT <- fortify(F252_mcp.EM, region = "id")
F252_mcp.DRYT <- fortify(F252_mcp.DRY, region = "id")
F252_mcp.MONT <- fortify(F252_mcp.MON, region = "id")
F252_mcp.PMT <- fortify(F252_mcp.PM, region = "id")

F36_mcp.EMT <- fortify(F36_mcp.EM, region = "id")
F36_mcp.DRYT <- fortify(F36_mcp.DRY, region = "id")
F36_mcp.MONT <- fortify(F36_mcp.MON, region = "id")
F36_mcp.PMT <- fortify(F36_mcp.PM, region = "id")

F66_mcp.EMT <- fortify(F66_mcp.EM, region = "id")
F66_mcp.DRYT <- fortify(F66_mcp.DRY, region = "id")
F66_mcp.MONT <- fortify(F66_mcp.MON, region = "id")
F66_mcp.PMT <- fortify(F66_mcp.PM, region = "id")

mcp.shift.TEST5 <- ggplot() +
  geom_polygon(data=F114_mcp.EMT, aes(x=F114_mcp.EMT$long, y=F114_mcp.EMT$lat),
               alpha=0.1,colour="blue",linetype=2) +
  geom_polygon(data=F114_mcp.DRYT, aes(x=F114_mcp.DRYT$long, y=F114_mcp.DRYT$lat),
               alpha=0.1,colour="red",linetype=3) +
  geom_polygon(data=F114_mcp.MONT, aes(x=F114_mcp.MONT$long, y=F114_mcp.MONT$lat),
               alpha=0.1,colour="green",linetype=4) +
  geom_polygon(data=F114_mcp.PMT, aes(x=F114_mcp.PMT$long, y=F114_mcp.PMT$lat),
               alpha=0.1,colour="black",linetype=5) +
  geom_polygon(data=F137_mcp.EMT, aes(x=F137_mcp.EMT$long, y=F137_mcp.EMT$lat),
               alpha=0.1,colour="blue",linetype=2) +
  geom_polygon(data=F137_mcp.DRYT, aes(x=F137_mcp.DRYT$long, y=F137_mcp.DRYT$lat),
               alpha=0.1,colour="red",linetype=3) +
  geom_polygon(data=F137_mcp.MONT, aes(x=F137_mcp.MONT$long, y=F137_mcp.MONT$lat),
               alpha=0.1,colour="green",linetype=4) +
  geom_polygon(data=F137_mcp.PMT, aes(x=F137_mcp.PMT$long, y=F137_mcp.PMT$lat),
               alpha=0.1,colour="black",linetype=5) +
  geom_polygon(data=F147_mcp.EMT, aes(x=F147_mcp.EMT$long, y=F147_mcp.EMT$lat),
               alpha=0.1,colour="blue",linetype=2) +
  geom_polygon(data=F147_mcp.DRYT, aes(x=F147_mcp.DRYT$long, y=F147_mcp.DRYT$lat),
               alpha=0.1,colour="red",linetype=3) +
  geom_polygon(data=F147_mcp.MONT, aes(x=F147_mcp.MONT$long, y=F147_mcp.MONT$lat),
               alpha=0.1,colour="green",linetype=4) +
  geom_polygon(data=F147_mcp.PMT, aes(x=F147_mcp.PMT$long, y=F147_mcp.PMT$lat),
               alpha=0.1,colour="black",linetype=5) +
  # geom_polygon(data=F252_mcp.EMT, aes(x=F252_mcp.EMT$long, y=F252_mcp.EMT$lat),
  #              alpha=0.1,colour="black",linetype=2) +
  # geom_polygon(data=F252_mcp.DRYT, aes(x=F252_mcp.DRYT$long, y=F252_mcp.DRYT$lat),
  #              alpha=0.1,colour="black",linetype=3) +
  # geom_polygon(data=F252_mcp.MONT, aes(x=F252_mcp.MONT$long, y=F252_mcp.MONT$lat),
  #              alpha=0.1,colour="black",linetype=4) +
  # geom_polygon(data=F252_mcp.PMT, aes(x=F252_mcp.PMT$long, y=F252_mcp.PMT$lat),
  #              alpha=0.1,colour="black",linetype=5) +
  geom_polygon(data=F36_mcp.EMT, aes(x=F36_mcp.EMT$long, y=F36_mcp.EMT$lat),
               alpha=0.1,colour="blue",linetype=2) +
  geom_polygon(data=F36_mcp.DRYT, aes(x=F36_mcp.DRYT$long, y=F36_mcp.DRYT$lat),
               alpha=0.1,colour="red",linetype=3) +
  geom_polygon(data=F36_mcp.MONT, aes(x=F36_mcp.MONT$long, y=F36_mcp.MONT$lat),
               alpha=0.1,colour="green",linetype=4) +
  geom_polygon(data=F36_mcp.PMT, aes(x=F36_mcp.PMT$long, y=F36_mcp.PMT$lat),
               alpha=0.1,colour="black",linetype=5) +
  geom_polygon(data=F66_mcp.EMT, aes(x=F66_mcp.EMT$long, y=F66_mcp.EMT$lat),
               alpha=0.1,colour="blue",linetype=2) +
  geom_polygon(data=F66_mcp.DRYT, aes(x=F66_mcp.DRYT$long, y=F66_mcp.DRYT$lat),
               alpha=0.1,colour="red",linetype=3) +
  geom_polygon(data=F66_mcp.MONT, aes(x=F66_mcp.MONT$long, y=F66_mcp.MONT$lat),
               alpha=0.1,colour="green",linetype=4) +
  geom_polygon(data=F66_mcp.PMT, aes(x=F66_mcp.PMT$long, y=F66_mcp.PMT$lat),
               alpha=0.1,colour="black",linetype=5) +
  theme_bw() +
  labs(x="Easting (m)", y="Northing (m)") +
  labs(caption = "Figure 6 |  SC seasonal home range shifts of five lizards. All seasonal polygons stay relatively stable with \n considerable overlap and without any major shifts.")+
  theme(plot.caption = element_text(hjust = 0,lineheight = 0.9))+
  theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5))

mcp.shift.TEST5

TABLE group means of seasonal home ranges between the two populations averaged across sex

seasonal<-read.csv("SC_Seasonal_Data.csv")

library(Rmisc)

SEAS_GRP_Means <- summarySE(seasonal, measurevar="Home_Range_100mcp",
                            groupvars=c("Environment","Season"), na.rm = TRUE)

# SEAS_GRP_Means
kable(SEAS_GRP_Means, format = "pandoc", caption = 'Table 6 | Group means of seasonal home ranges between Stone Canyon (subsidized) and Owl Head Buttes (non-subsidized). These means are averaged across sex.')

RM-ANOVA for seasonal home ranges between environments

library(lme4)
library(readr)
library(lmerTest)
# seasonal<-read.csv("SC_Seasonal_Data.csv")

RM.mod.Season <- lmer(Home_Range_100mcp~Environment+Season+Sex+N+Environment*Season+(1|Gila), 
                      data=seasonal)
summary(RM.mod.Season)

ANOVA table of seasonal HRs between envs.

anova(RM.mod.Season)

TABLE of seasonal home ranges by sex between the two populations

SEAS_GRP_TEST <- summarySE(seasonal, measurevar="Home_Range_100mcp",
                           groupvars=c("Environment","Season","Sex"), na.rm = TRUE)

# SEAS_GRP_Means
kable(SEAS_GRP_TEST, format = "pandoc", caption = 'Table 7 | Seasonal home range means between Stone Canyon (subsidized) and Owl Head Buttes (non-subsidized) popuations for males and females. These are raw means before being adjusted for environment, season, sex, and sample size.')

figures for raw seasonal home ranges between the two populations

Figures Adjusted EMMs of seasonal home range between the two populations

Collective grid of raw and adjusted seasonal home ranges

ggarrange(raw.seasonal, adj.seasonal, labels = c("A", "B"),
          nrow = 2)

Post hoc analyses of seasonal home ranges

Pairwise of each season between populations, overaged over levels of sex

emm_s.t <- emmeans(RM.mod.Season, pairwise ~ Environment | Season)
emm_s.t

Graphical comparisons

plot(emm_s.t, comparisons = TRUE)

Pairwise between seasons within each popultion

emm_s.t4 <- emmeans(RM.mod.Season, pairwise ~ Season | Environment)
emm_s.t4

Graphical Comps

plot(emm_s.t4, comparisons = TRUE)

Pairwise between sexes of each season of the subsidized population

sub <- subset(seasonal, Environment == "subsidized")

RM.mod.Sub <- lmer(Home_Range_100mcp~Season+Sex+N+Season*Sex+(1|Gila), data=sub)

emm_s.t5 <- emmeans(RM.mod.Sub, pairwise ~ Sex | Season)
emm_s.t5 

Graphical Comps

plot(emm_s.t5, comparisons = TRUE)

Pairwise between sexes of each season of the non-subsidized population

nonsub <- subset(seasonal, Environment == "nonsubsidized")
View(nonsub)
RM.mod.NSub <- lmer(Home_Range_100mcp~Season+Sex+N+Season*Sex+(1|Gila), data=nonsub)

emm_s.t6 <- emmeans(RM.mod.NSub, pairwise ~ Sex | Season)
emm_s.t6 

Graphical Comps

plot(emm_s.t6, comparisons = TRUE)